/// Script to produce relative risk graph and estimates using nlcom command and coefplot package
/// Fabien Cottier
/// SPEI DATA


///////// Graph immediate SPEI effect only(SPEI Y0)

* write programe for relative risk estimation // margins cannot produce rr estimates for cont. variables
program define rr_est	
	local models ""
	local x=-.5
	
	* Loop to create multiple nlcom statement
	forvalues j=1(1)11 {

		if "`2'"=="linear"{
			local models `"`models' (_b[wmean_spei]*`x')"'
		}
		
		else if "`2'"=="quadratic"{
			local models `"`models' (_b[wmean_spei]*`x'+_b[wmean_spei#wmean_spei]*(`x')^2)"'
		}
		
		local x=`x'+1/10
	}

	* estimates relative risk using nlcom
	est resto `1' 
	local df_est=e(df_r)
	nlcom `models', post df(`df_est') 
	estimates store est_rr
	
	* produce plot using coefplot
	coefplot (est_rr), ///
		eform yscale(range(.75 1.75)) ylabel(.75 "-25%" 1 "0" 1.25 "+25%" 1.5 "+ 50%" 1.75 "+75%") xlabel(1 "-0.5" 3.5 "-0.25" 6 "0" 8.5 "0.25" 11 "0.5") ///
		ytitle(Relative change in irr. migration (%)) xtitle(Standardized Precipitation Evapotranspiration Index) vertical recast(line) ciopts(recast(rarea)) yline(1) 
end




///////// Graph effects for multiple lags SPEI (SPEI Y0, Y-1 Y-2)


* write programe for relative risk estimation // margins cannot produce rr estimates for cont. variables
program define rr_est_MuL	// MuL for multiple lag
	local models_L0 ""  
	local models_L1 ""  
	local models_L2 ""  
	local x=-.5
	* Loop to create multiple nlcom statement
	forvalues j=1(1)11 {
	
		if "`2'"=="linear"{
			local models_L0 `"`models_L0' (_b[wmean_spei]*`x')"'
			local models_L1 `"`models_L1' (_b[L4.wmean_spei]*`x')"'
			local models_L2 `"`models_L2' (_b[L8.wmean_spei]*`x')"'
		}
	
		else if "`2'"=="quadratic"{
			local models_L0 `"`models_L0' (_b[wmean_spei]*`x'+_b[wmean_spei#wmean_spei]*(`x')^2)"'
			local models_L1 `"`models_L1' (_b[L4.wmean_spei]*`x'+_b[L4.wmean_spei#L4.wmean_spei]*(`x')^2)"'
			local models_L2 `"`models_L2' (_b[L8.wmean_spei]*`x'+_b[L8.wmean_spei#L8.wmean_spei]*(`x')^2)"'
		}
		
		local x=`x'+1/10
	}
	
	* estimates relative risk using nlcom
	est resto `1'
	local df_est=e(df_r)
	nlcom `models_L0', post df(`df_est')
	estimates store est_rr_L0
	est resto `1'
	nlcom `models_L1', post df(`df_est')
	estimates store est_rr_L1
	est resto `1'
	nlcom `models_L2', post df(`df_est')
	estimates store est_rr_L2
	
	* graph
	coefplot (est_rr_L0), bylabel(Immediate effect) || ///
		(est_rr_L1), bylabel(Lag effect (Year-1)) || ///
		(est_rr_L2), bylabel(Lag effect (Year-2)) ||, ///
		eform yscale(range(.75 1.75)) ylabel(.75 "-25%" 1 "0" 1.25 "+25%" 1.5 "+ 50%" 1.75 "+75%") xlabel(1 "-0.5" 3.5 "-0.25" 6 "0" 8.5 "0.25" 11 "0.5") ///
		ytitle(Relative change in irr. migration (%)) xtitle(Standardized Precipitation Evapotranspiration Index) vertical recast(line) ciopts(recast(rarea)) yline(1)
		
end



///////// Split sample graph effects SPEI (SPEI Y0)


* write programe for relative risk estimation // margins cannot produce rr estimates for cont. variables
program define rr_est_Sp	// MuL_Sp for multiple lag with split sample
	local models_L0 ""  
	local x=-.5
	* Loop to create multiple nlcom statement
	forvalues j=1(1)11 {
	
		if "`3'"=="linear"{
			local models_L0 `"`models_L0' (_b[wmean_spei]*`x')"'
		}
	
		else if "`3'"=="quadratic"{
			local models_L0 `"`models_L0' (_b[wmean_spei]*`x'+_b[wmean_spei#wmean_spei]*(`x')^2)"'
		}
		
		local x=`x'+1/10
	}
	
	* estimates relative risk using nlcom // With High labor agriculture
	est resto `1'
	local df_est=e(df_r)
	nlcom `models_L0', post df(`df_est')
	estimates store est_rr_L0H

		
	* estimates relative risk using nlcom // With Low labor agriculture
	est resto `2'
	local df_est2=e(df_r)
	nlcom `models_L0', post df(`df_est2')
	estimates store est_rr_L0L
	
	* graph
	coefplot (est_rr_L0H), bylabel(High Agriculture - Immediate effect) || ///
		(est_rr_L0L), bylabel(Low Agriculture - Immediate effect) ||, ///
		eform yscale(range(.75 1.75)) ylabel(.75 "-25%" 1 "0" 1.25 "+25%" 1.5 "+ 50%" 1.75 "+75%") xlabel(1 "-0.5" 3.5 "-0.25" 6 "0" 8.5 "0.25" 11 "0.5") ///
		ytitle(Relative change in irr. migration (%)) xtitle(Standardized Precipitation Evapotranspiration Index) vertical recast(line) ///
		ciopts(recast(rarea)) yline(1)  byopts(row(3))
		
end




 
